
time_array_min = (datenum(1995,1,1,0,0:9468000-1,30));
a = length(time_array_min);
nc_fname = '/Users/manojnair/Downloads/DED_1995_2012.nc'
[x_data, y_data, z_data , X_ID, Y_ID, Z_ID, obj, ncid] = read_geomag_netcdf(nc_fname, 0, a, 0);

%fit spline
    
    L = isnan(x_data);
    data_array = x_data(~L);
    time_array = time_array_min(~L);
    b1 = min(time_array):365:max(time_array)+10;
    if length(b1) > 1,
        sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)));
        v=ppval(time_array_min,sp);
    else,%data length <=1 year
        sp=robustfit(time_array,data_array);
        v= sp(1) + sp(2) * time_array_min;
    end;
    subplot(311);
    plot(time_array_min,x_data);
    hold on;
    plot(time_array_min(~L),v(~L),'r');
    ylabel('X');
    axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
    datetick('x','keeplimits');
    title('DED');
    % write sec removed z values
    x_data = x_data - v';
    
    L = isnan(y_data);
    data_array = y_data(~L);
    time_array = time_array_min(~L);
    b1 = min(time_array):365:max(time_array)+10;
    if length(b1) > 1,
        sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)));
        v=ppval(time_array_min,sp);
    else,%data length <=1 year
        sp=robustfit(time_array,data_array);
        v= sp(1) + sp(2) * time_array_min;
    end;
    subplot(312);
    plot(time_array_min,y_data);
    hold on;
    plot(time_array_min(~L),v(~L),'r');
    ylabel('Y');
    axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
    datetick('x','keeplimits');
    
    % write sec removed z values
    y_data = y_data - v';
    
    
    L = isnan(z_data);
    data_array = z_data(~L);
    time_array = time_array_min(~L);
    b1 = min(time_array):365:max(time_array)+10;
    if length(b1) > 1,
        sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)));
        v=ppval(time_array_min,sp);
    else,%data length <=1 year
        sp=robustfit(time_array,data_array);
        v= sp(1) + sp(2) * time_array_min;
    end;
    subplot(313);
    plot(time_array_min,z_data);
    hold on;
    plot(time_array_min(~L),v(~L),'r');
    ylabel('Z');
    axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
    datetick('x','keeplimits');
    
    % write sec removed z values
    z_data = z_data - v';
    
    
    % Now store this data back into the same file. The idea is that the
    % file is then renamed into Secular_Removed.nc
    
ncid = netcdf.open(nc_fname,'NC_WRITE');
x_data(isnan(x_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
y_data(isnan(y_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
z_data(isnan(z_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)

%multiplied by 10 to preserve the decimal points
netcdf.putVar(ncid, X_ID, 0, a,  int32 (x_data .* 10) );
netcdf.putVar(ncid, Y_ID, 0, a,  int32 (y_data .* 10) );
netcdf.putVar(ncid, Z_ID, 0, a,  int32 (z_data .* 10) );
netcdf.close(ncid);

    
    
    
